Image phenotyping platforms

Image phenotyping platforms are a fast and non-invasive tool for phenotyping plant under lab and field conditions. The most common are sensors for imaging. The reflected light captured in these images can be used to draw inference about many traits, including:

• Geometric traits (i.e. plant height, leaf area index, lodging, crop canopy cover)
• Canopy spectral texture (spectral features)
• Physiological traits (i.e., chlorophyll, biomass, pigment content, photosynthesis)
• Abiotic/biotic stress indicators (i.e., stomatal conductance, canopy temperature difference, leaf water potential, senescence index)
• Nutrients (nitrogen concentration, protein content)
• Yield

Image segmentation

Image segmentation is the process of partitioning a digital image into multiple segments (sets of pixels or image objects).

Thresholding

Thresholding is the simplest method of image segmentation. This method is based on a clip-level (or a threshold value) to turn a gray-scale image into a binary image (Mask).

FIELDimageR

FIELDimageR is a R package to analyze images from research fields, and allows to:

• Crop the image
• Remove the background
• Build vegetation indices
• Rotate the image
• Build the plot shapefile
• Extract information for each plot

1) Example 01

Counting green seeds in Soybean

Alexandre Abdão(Embrapa Maize & Sorghum)

• 4 pictures: soybean samples with different green seeds percentage

• Number of green seeds (0, 1, 6, and 11)

• Donwload: soybean.zip

Required packages

• FIELDimageR

# Installing 
install.packages("sp")
install.packages("raster")
install.packages("rgdal")
install.packages("ggplot2")
install.packages("agricolae")
install.packages("reshape2")
install.packages("devtools")
install.packages("parallel")
install.packages("foreach")
install.packages("doParallel")
install.packages("lme4")
devtools::install_github("filipematias23/FIELDimageR")

# Necessary packages 
library(FIELDimageR)
library(raster)
library(ggplot2)
library(agricolae)
library(reshape2)
library(lme4)

• Uploading one image as example and decreasing the resolution

EX1<-stack("soybean/11.jpg")
EX1<-aggregate(EX1, fact= 4)
plotRGB(EX1)

# extent=T (The whole image area will be the shapefile)
EX.shapeFile<-fieldPolygon(EX1,extent = T)

# Select one index to identify leaves and remove the background
EX1.I1<- fieldIndex(mosaic = EX1,index = c("SI","BGI","BI"))

# Thresholding
dev.off()
par(mfrow=c(1,2))
hist(EX1.I1$BGI)
plot(EX1.I1$BGI)

# Removing the background
EX1.R1<- fieldMask(mosaic = EX1, index = "BGI",
                   cropValue = 0.7,
                   cropAbove = T)

# Counting the total number of seeds
EX.P.Total<-fieldCount(mosaic = EX1.R1$mask, 
                       fieldShape = EX.shapeFile$fieldShape, 
                       minSize = 0.1,
                       cex = 1.5,
                       na.rm = T)

# Select one index to identify green seeds
EX1.I2<- fieldIndex(mosaic = EX1.R1$newMosaic,index = c("SI","BGI","BI"))

#BI index
plot(EX1.I2$BI)

# Selecting green seeds
EX1.R2<- fieldMask(mosaic = EX1.R1$newMosaic, 
                   index = "BI",
                   #myIndex = "Blue",
                   cropValue = 130,
                   cropAbove = T)

# Counting the number of green seeds
EX.Green<-fieldCount(mosaic = EX1.R2$mask, 
                     fieldShape = EX.shapeFile$fieldShape, 
                     minSize = 0.07,
                     cex = 1.5,
                     na.rm = T)

# Joying information
data.frame(Total=EX.P.Total$fieldCount,
           Green=EX.Green$fieldCount,
           Percentage=round(EX.Green$fieldCount/EX.P.Total$fieldCount,2))

################
### Parallel ###
################

# Required packages
library(parallel)
library(foreach)
library(doParallel)

# Images names (folder directory: "./soybean/")
pics<-list.files("./soybean/")

# Number of cores
n.core<-2

# Starting parallel
cl <- makeCluster(n.core, output = "")
registerDoParallel(cl)
EX.Table.Parallel <- foreach(i = 1:length(pics), .packages = c("raster","FIELDimageR"), 
                               .combine = rbind) %dopar% {
EX1<-stack(paste("./soybean/",pics[i],sep = ""))
EX1<-aggregate(EX1, fact= 4)
EX.shapeFile<-fieldPolygon(EX1,extent = T,
                           plot = F)
EX1.R1<- fieldMask(mosaic = EX1, index = "BGI",
                   cropValue = 0.7,
                   cropAbove = T,
                   plot = F)
EX.P.Total<-fieldCount(mosaic = EX1.R1$mask, 
                       fieldShape = EX.shapeFile$fieldShape, 
                       minSize = 0.1,
                       cex = 1.5,
                       na.rm = T)
EX1.R2<- fieldMask(mosaic = EX1.R1$newMosaic, 
                   index = "BI",
                   cropValue = 130,
                   cropAbove = T,
                   plot = F)
EX.Green<-fieldCount(mosaic = EX1.R2$mask, 
                     fieldShape = EX.shapeFile$fieldShape, 
                     minSize = 0.07,
                     cex = 1.5,
                     na.rm = T)
data.frame(Total=EX.P.Total$fieldCount,
           Green=EX.Green$fieldCount,
           Percentage=round(EX.Green$fieldCount/EX.P.Total$fieldCount,2))
}
rownames(EX.Table.Parallel)<-pics
EX.Table.Parallel 

2) Example 02

Remote sensing (Potato Breeding)

Jeffrey Endelman and Filipe Matias (UW-Madison)

• Donwload: odm_orthophoto.tif

# Uploading file (odm_orthophoto.tif):
EX1<-stack("odm_orthophoto.tif")
plotRGB(EX1)

# Rotating the image with theta 2.3:
EX1.Rotated<-fieldRotate(mosaic = EX1[[-4]],theta = 2.3)

# Removing soil (example using HUE):
EX1.RemSoil<- fieldMask(mosaic = EX1.Rotated)

• Evaluating each trial per time:

# Trial: 01
EX1.Crop.1 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 

• Dataset - Field notes - phenotype (Download: ‘EX1_Data.xlsx’)

# Dataset - Field notes - phenotype (.exel or .txt):
Data1<-read_excel('EX1_Data.xlsx',1)

# Field map ID identification (ID="Plot"):
Map1<-fieldMap(fieldPlot = Data1$Plot,fieldColumn = Data1$Column, fieldRow = Data1$Row,decreasing = T)
Map1

# Building the plot shapefile (ncols = 14 and nrows = 10)
x11()
EX1.Shape.1<-fieldShape(mosaic = EX1.Crop.1, ncols = 14, nrows = 10, fieldData = Data1, ID = "Plot", fieldMap = Map1)

# Trial: 02
EX1.Crop.2 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data2<-read_excel('EX1_Data.xlsx',2)
Map2<-fieldMap(fieldPlot = Data2$Plot,fieldColumn = Data2$Column, fieldRow = Data2$Row,decreasing = T)
EX1.Shape.2<-fieldShape(mosaic = EX1.Crop.2, ncols = 14, nrows = 9, fieldData = Data2, ID = "Plot", fieldMap = Map2)

# Trial: 03
EX1.Crop.3 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data3<-read_excel('EX1_Data.xlsx',3)
Map3<-fieldMap(fieldPlot = Data3$Plot,fieldColumn = Data3$Column, fieldRow = Data3$Row,decreasing = T)
EX1.Shape.3<-fieldShape(mosaic = EX1.Crop.3, ncols = 14, nrows = 13, fieldData = Data3, ID = "Plot", fieldMap = Map3)

# Trial: 04
EX1.Crop.4 <- fieldCrop(mosaic = EX1.RemSoil$newMosaic) 
Data4<-read_excel('EX1_Data.xlsx',4)
Map4<-fieldMap(fieldPlot = Data4$Plot,fieldColumn = Data4$Column, fieldRow = Data4$Row,decreasing = T)
EX1.Shape.4<-fieldShape(mosaic = EX1.Crop.4, ncols = 14, nrows = 10, fieldData = Data4, ID = "Plot", fieldMap = Map4)

# Combining shapefiles
EX1.Shape<-rbind(EX1.Shape.1$fieldShape,EX1.Shape.2$fieldShape,EX1.Shape.3$fieldShape,EX1.Shape.4$fieldShape)
plot(EX1.Shape)

# Building indices ("NGRDI","BGI","GLI","SCI")
EX1.Indices<- fieldIndex(mosaic = EX1.RemSoil$newMosaic,
                         index = c("NGRDI","BGI"), 
                         myIndex = c("(Red-Blue)/Green"))

plot(EX1.Indices$myIndex)
plot(EX1.Shape,add=T)

# Extracting Data
EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI","myIndex")]],
                     fieldShape = EX1.Shape,
                     buffer = -0.05,
                     n.core = 3)
NewData<-EX1.Info$fieldShape@data
NewData

# Making map plots
fieldPlot(fieldShape=EX1.Info$fieldShape,
          fieldAttribute="NGRDI", 
          mosaic=EX1.Indices, 
          color=c("red","blue"), 
          alpha = 0.4,
          round = 2)

• Estimation of plant height using a high throughput phenotyping platform. (Donwload: dsm.tif)

# Uploading files from vegetative growth (dsm_1.tif):
DSM <- stack("dsm.tif")
plot(DSM)

# Rotating the image using the same theta (2.3):
DSM.R<-fieldRotate(DSM, theta = 2.3)

# Removing the soil using mask from step 4:
DSM.S <- fieldMask(DSM.R, mask = EX1.RemSoil$mask)

# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(DSM.S$newMosaic, fieldShape = EX1.Info$fieldShape, fun = "mean",n.core = 3)
EPH$plotValue

# Correlation
NewData2<-EPH$fieldShape@data[,c("layer","NGRDI","BGI","myIndex")]
cor1<-correlation(NewData2)
cor1$correlation
cor1$pvalue

# Regression
NewData3<-EPH$fieldShape@data[,c("Trial","Name","Block","Row","Column","layer","NGRDI")]
NewData3$Trial<-as.factor(NewData3$Trial)
NewData3$Name<-as.factor(NewData3$Name)
NewData3$Block<-as.factor(NewData3$Block)
NewData3$Row<-as.factor(NewData3$Row)
NewData3$Column<-as.factor(NewData3$Column)
NewData3$layer<-as.numeric(as.character(NewData3$layer))
NewData3$NGRDI<-as.numeric(as.character(NewData3$NGRDI))

ggplot(NewData3,aes(x=layer,y=NGRDI, col=Trial))+
  facet_wrap(~Trial, scales = "fixed",ncol = 1)+
  geom_point() +
  geom_smooth(method=lm)+
  theme_bw()

fieldPlot(fieldShape=EPH$fieldShape,
          fieldAttribute="layer", 
          mosaic=EX1.Indices, 
          color=c("red","blue"), 
          alpha = 0.4,
          round = 2)

• Heritability

# NGRDI

str(NewData3)

mod.T1<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)

mod.T2<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)

mod.T3<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)

mod.T4<-lmer(NGRDI~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)

NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
                     H2=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))

ggplot(NewData4,aes(x=Trail,y=H2,fill=Trail))+
  geom_bar(stat="identity")+
  geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
  theme_bw()

# Estimated plant height 

mod.T1<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T1",])
H2.T1<-as.data.frame(VarCorr(mod.T1))$vcov[1]/sum(as.data.frame(VarCorr(mod.T1))$vcov)

mod.T2<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T2",])
H2.T2<-as.data.frame(VarCorr(mod.T2))$vcov[1]/sum(as.data.frame(VarCorr(mod.T2))$vcov)

mod.T3<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T3",])
H2.T3<-as.data.frame(VarCorr(mod.T3))$vcov[1]/sum(as.data.frame(VarCorr(mod.T3))$vcov)

mod.T4<-lmer(layer~Row+Column+(1|Name),NewData3[NewData3$Trial=="T4",])
H2.T4<-as.data.frame(VarCorr(mod.T4))$vcov[1]/sum(as.data.frame(VarCorr(mod.T4))$vcov)

NewData4<-data.frame(Trail=factor(c("T1","T2","T3","T4")),
                     EPH=as.numeric(c(H2.T1,H2.T2,H2.T3,H2.T4)))

ggplot(NewData4,aes(x=Trail,y=EPH,fill=Trail))+
  geom_bar(stat="identity")+
  geom_text(aes(label=round(EPH,2)), vjust=1.6, color="white", size=6)+
  theme_bw()